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A computational investigation was performed to assess the effectiveness of Detached- 
Eddy Simulation (DES) as a tool for predicting icing effects. The AVUS code, a public 
domain flow solver, was employed to compute solutions for an iced wing configuration using 
DES and steady Reynolds Averaged Navier-Stokes (RANS) equation methodologies. The 
wing section considered here was a business jet airfoil (GLC305) with a 22.5-minute glaze 
ice accretion (944-ice shape). The section was extruded to form a rectangular planform. 
The model was mounted between two walls so no tip effects were considered. The numerical 
results were validated by comparison with experimental data for the same configuration. 
The time-averaged DES computations showed some improvement in lift and drag results 
near stall when compared to steady RANS results. However, comparisons of the flow field 
details did not show the level of agreement suggested by the integrated quantities. More 
validation is needed to determine what role DES can play as part of an overall icing effects 
prediction strategy. 


I. Introduction 

It is well known and appreciated that the accretion of ice on lifting surfaces can significantly degrade the 
performance and handling characteristics of affected aircraft. 1 Less well understood, however, are the details 
of the complex flow fields associated with these configurations. Depending on the ice shape, these flow fields 
may exhibit extensive regions of unsteady separated flow. Typically, ice accretions on straight wings are 
categorized as “rime,” “horn,” or “spanwise ridge,” which are more or less two-dimensional in nature, or 
“roughness,” which is the most three-dimensional of the shapes. 1 The flow field for a horn ice accretion, 
which is the subject of this effort, is dominated by a separation bubble downstream of the horn. The 
flow separates because it is unable to recover the necessary pressure in the boundary layer after it expands 
around the horn. The ensuing free shear layer rolls up to form vortices that are convected downstream. The 
development of this unsteady flow field is critical to the resulting flow near maximum lift. The separation 
bubble has a large global effect and is similar to the long bubble defined by Tani. 2 

Many researchers have turned to computational fluid dynamics (CFD) simulations to investigate iced 
wing flow fields. Beginning with Potapczuk, 3 numerous results from icing effects studies using the Reynolds 
averaged Navier-Stokes (RANS) equations have been reported in the literature. Chung, et at analyzed the 
flow around the ice contaminated wing surfaces of a turbo-prop aircraft. Their two- and three-dimensional 
analyses of a wing with a spanwise ridge accretion employed a steady RANS approach and were performed 
to provide insight into the aerodynamics that may have led to a loss of control of the aircraft. Reported 
discrepancies between the two- and three-dimensional results were attributed to a lack of grid resolution 
for the three-dimensional cases considered. Dunn and Loth 5 studied the effects of simulated spanwise ice 
shapes on the aerodynamic performance of two-dimensional wing sections using a RANS solver developed for 
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unstructured meshes. A forward-facing quarter round was employed as the simulated ice shape. Additionally, 
they employed a solution adaptive mesh in which the mesh was refined in the separated shear layer. Their 
results showed that, even with adaptive mesh refinement, the pressure recovery was not correctly predicted. 
This discrepancy was attributed to the inability of the steady RANS simulation to model the experimentally 
observed unsteadiness that occurs in flow fields of this type. Pan, Loth, and Bragg 6 performed RANS 
simulations for flow about airfoils with ridge ice shapes and leading edge ice shapes and compared the 
results to experimental data. Favorable comparisons between predicted force data and experiments were 
obtained up to, but not including, the stall condition, which is dominated by unsteady flow. 

Recently, researchers have begun employing unsteady simulations to investigate separated flow fields. 
Candidate approaches for computing such unsteady flow fields include unsteady RANS (URANS), direct 
numerical simulation (DNS), large-eddy simulation, (LES), and detached eddy simulation (DES). Spalart 7 
contends that URANS is both “ambiguous and flawed” and that “its quantitative performance can be quite 
poor.” DNS, while attractive from the standpoint that no modeling is required, and LES are currently 
impractical for realistic configurations due to their computational expense. DES, however, appears to have 
the potential to represent a viable near-term approach. DES is a hybrid RANS/LES approach in which a 
RANS turbulence model is employed in attached thin shear layers and an LES-type model is employed in 
regions away from the wall. This approach exploits the ability of RANS simulations to efficiently model 
high Reynolds number attached boundary layer flows and the ability of LES to model geometry depen- 
dent, unsteady three-dimensional flows. 8 DES has been employed in numerous simulations of complex flow 
fields. 9 - 10 ’ 11 DES has also been applied to the icing effects problem. Kumar and Loth 12 employed a DES 
technique to predict the unsteady flow around an NLF0414 airfoil section with a synthetic ice shape located 
on the upper surface at 3.4% chord. They also reported on computations for a rectangular three-dimensional 
wing using the same section. Their three-dimensional results showed some improvement relative to the 
steady RANS computations. Pan and Loth 13 performed a DES for an NACA23012 airfoil with forward- 
facing quarter-round simulated ice accretion. Again, some improvement between the predicted results and 
experimental data was reported. 

We previously reported results from steady RANS simulations for a horn ice shape/wing configuration. 14 
The basic ice shape considered was a 22.5-minute glaze ice accretion on a GLC305 airfoil, which is denoted 
as the two-dimensional 944- ice shape. 15 The airfoil section/ice shape was extruded to form a rectangular 
planform wing with an aspect ratio of unity. The flow fields associated with this configuration have been 
studied experimentally and the results reported in the literature. 16 ’ 17 These data include force and pressure 
measurements as well as mean velocity measurements and RMS fluctuations of the streamwise and transverse 
velocity components obtained using a split-hot-film probe. The model was mounted between the tunnel walls, 
so no tip effects were considered. We compared steady RANS results computed using the non-commercial 
version of the Cobaltgo code 18 (now called AVUS) with the Spalart-Allmaras turbulence model 19 with the 
experimental data. In general, the steady RANS computations underpredicted both the lift and drag as 
the angle of attack was increased. Additionally, the extent of the separated region was overpredicted even 
at lower angles of attack. As noted above for other RANS computations, the pressure recovery was not in 
agreement with the experimental data. In a related effort, Chi, et al . 20 reported on the results of a series 
of two-dimensional steady RANS simulations for the GLC305/944-ice shape, as well as a rime ice shape, 
using a variety of turbulence models. The two-dimensional results for the 944-ice shape were qualitatively 
very similar to the three-dimensional RANS computations 14 regardless of the turbulence model employed. 
Additionally, relatively fine meshes were employed to discretize the computational domain. Therefore, we 
believe the discrepancies reported by Chi, et al. 20 were not related to mesh effects. 

In this paper we present the results of steady RANS computations and unsteady DES computations 
performed for the rectangular planform, extruded GLC3Q5/944-ice shape wing. We first provide a brief 
discussion of the DES approach. Numerical considerations regarding how the mesh spacing and time step size 
are selected along with the mesh generation and solution procedures are discussed next. We present sample 
results from the two- and three-dimensional computations and compare these results with experimental data. 


II. Problem Definition 

The specific airfoil section considered during this effort was the GLC305 airfoil section with the 22.5- 
minute glaze ice accretion which is denoted as the two-dimensional 944-ice shape. 15 The ice shape was 
extruded to form a rectangular planform wing with an aspect ratio of unity to match the configuration used 
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in the test program 16 * 17 and is referred to here as the “extruded wing.” The airfoil model was supported 
horizontally across the width of the test section between two circular end-plates. These end-plates were 
flush mounted with the sidewalls and rotated to provide the angle of attack adjustment. Each end-plate 
was equipped with a porous section for sidewall boundary layer control. Numerical solutions were compared 
with experimental data for the following conditions: M=0.12, Re/L=3.8xl0 6 /m which, with a chord length 
of 0.9144m, yields Re=3.5xl0 6 and angles of attack from 0 deg to 6 deg in 2 deg increments. The extruded 
944-ice shape cases considered here correspond to Run 41 in the experimental data. 16 * 17 

III. Detached Eddy Simulation 

DES was originally formulated as a modification to the Spalart-Allmaras turbulence model. 19 ’ 8 The 
basic idea is to employ a RANS-type turbulence model for thin, attached shear layers and an LES model 
in separated regions. This approach exploits the ability of RANS simulations to efficiently model attached 
boundary layer flows and the ability of LES to model unsteady, geometry-dependent, three-dimensional flow 
fields. RANS-type models perform well for attached flows and less well for separated flows, while LES models 
perform well for separated flows but are prohibitively expensive to employ for computations involving thin 
boundary layers. We now provide a brief overview of the DES model as implemented in the Spalart-Allmaras 
turbulence model. More detailed discussions can be found in the literature. 8,10 ’ 11 ’ 21 

In the Spalart-Allmaras model, the wall destruction term is taken to be proportional to ( y/d ) 2 where 
£ is a working variable related to the turbulent viscosity and d is the distance to the nearest wall. If local 
equilibrium occurs and the production and destruction terms are in balance, the eddy viscosity becomes 
proportional to Sd 2 where S is the local strain rate. In the Smagorinsky LES model, 22 the subgrid-scale 
turbulent viscosity is proportional to the local strain rate and the local mesh spacing squared, i.e., SA? 
where A = max ( Ax , Ay, A z). Thus, if the distance to the wall d is replaced by the local mesh spacing A, 
the Spalart-Allmaras model will locally behave like a Smagorinsky LES model. To retain the RANS-type 
behavior in attached boundary layers, d is replaced by a new variable d = min(d, Cdes A) where Cdes 
is a constant. Note that when d « Cdes A, the model is in a RANS mode and models the average 
properties of attached flow turbulence. When d » Cdes A, the model is in LES mode and resolves eddies 
larger than some wavelength depending on the characteristics of the flow solved and the local mesh spacing. 
For unstructured meshes, the spacing A is taken to be the longest distance between the cell center and 
neighboring cell centers. 21 The implementation of the DES model has two immediate implications regarding 
the mesh: 

1. The mesh in the boundary layer should be highly anisotropic. In particular, the spacing along the 
surface should be larger than the local boundary layer thickness. This ensures that the model operates 
in a RANS mode near the boundary. 

2. There is no advantage to having an anisotropic mesh in the interior of the domain. The premise of 
LES is to filter out only eddies that are statistically isotropic 7 so that equal resolution in all directions 
is reasonable. 

Other implications regarding the discretization of the domain are discussed in the next section. 

IV. Numerical Considerations 

The prediction of the aerodynamic characteristics of an iced wing is a complex problem that involves 
several steps. We now describe each step in some detail. 

A. Geometry Modeling 

The GLC305/944-ice shape definition was used as input to ICEG2D, 23 a two-dimensional tool that automates 
geometry modeling and mesh generation for ice accretion predictions. ICEG2D redistributes points on the 
upper and lower surfaces of the defined airfoil using a curvature-based equidistribution algorithm to ensure 
a sufficient number of points are employed in regions of the airfoil surface with high curvature. A structured 
surface mesh was then extruded from the section definition. This surface mesh was then converted to a 
meshing-ready NURBS representation using the mesh generation software GUM-B. 24 
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Figure 1. Estimated DES zones for GLC305 airfoil section with two-dimensional 944-ice shape (extruded wing) 


B. Mesh Generation: Estimated DES Zones for an Iced Wing 

It is well known that the quality of the CFD solution depends on an appropriate mesh with sufficient mesh 
density in regions of high flow gradients. In an unsteady separated flow dominated by convecting vortices, 
the vortices must be resolved adequately along with the attached boundary layers. Indeed, as noted by 
Spalart, 7 “DES compounds the gridding difficulty by incorporating both types of turbulence treatment in the 
same field. ’ Spalart provides guidance as to how a domain should be discretized for DES modeling. He 
identifies several different regions with different meshing requirements: 

• An Euler region (ER) is a region that is free of turbulence and vorticity unless it is penetrated by a 
shockwave. The ER covers most of the domain and an isotropic mesh can be employed in this region. 

• A RANS region (RR) is primarily composed of the boundary layer where there is no LES content. 

- The viscous region (VR) is located within the RANS region and its meshing requirements are 
similar to those of standard RANS computations. 

— The RANS outer region (OR) should be discretized with a mesh in which the mesh spacing 
in the direction normal to the wall should not exceed once-tenth the thickness of the boundary 
layer. This is primarily an issue associated with numerical robustness. In practice, this rule is 
often violated. 

• An LES region (LR) contains vorticity and turbulence but is not a boundary layer. Significant LES 
content is present in an LR. 

- The focus region (FR) is the region near the body in which the separated turbulence must be 
well resolved. The mesh should be isotropic in the FR since the LES mode filters out eddies that 
are statistically isotropic. 

- The departure region (DR) is a transitional region between the FR and the ER. The DR does 
not need to be resolved as well as the FR. 

Here, regions are not distinguished by different equations being applied but by different priorities in the 
mesh spacing. An efficient mesh for any external flow can be designed with these concepts in mind, but not 
all are strict requirements. Figure 1 depicts the estimated DES zones discussed above as applied to the flow 
surrounding the 944-ice shape. 

C. Selection of Mesh Spacing in Focus Region 

The procedure described below follows that suggested by Spalart. 7 According to Spalart, a well-adjusted 
subgrid-scale model should allow energy cascade to the smallest eddies that can be resolved on the mesh. 
Therefore, for most CFD solvers, an eddy with a wavelength of A = 5Ao, where Aq is the local mesh spacing, 
will be active even though it cannot be highly accurate because it lacks the energy cascade to smaller eddies, 
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and is under the influence of eddy viscosity instead. In other words, if we want to adequately resolve an 
eddy with wavelength A, the mesh spacing should be A 0 = A/5. In the baseline mesh employed for DES 
computations, an eddy with a wavelength of 5% of the chord was selected, resulting in mesh spacing in the 
focus region of 1% of the chord. This wavelength was chosen because the height of the horn is approximately 
5% of the chord. Therefore, vortices with approximately the same spatial extent as the horn height can be 
resolved in the separated region downstream of the horn. As part of the comparison of the DES results with 
experimental data, the DES mesh was refined by a factor of two in the focus region. 

D. Selection of Time-Step Size 

Selecting the time-step size for DES computations is challenging. Spalart 19 suggests employing a local CFL 
number (based on the local flow velocity, the local mesh spacing, and the time step size) of unity in the 
focus regions (FR), that is -^^=1 where U 0 is the maximum flow velocity in the region. Computed RANS 
results indicate that the maximum flow velocity over the domain is approximately 40-50% greater than the 
freestream velocity. For the problem considered here with a freestream velocity of approximately 41.1m/s 
and a wing chord of 0.9144m, the computation for the time step size yields zU=0.15ms for the baseline 
mesh. For the refined mesh, ZU=0.075ms. However, as Spalart notes, 7 steps a factor of 3/2 or even 2 away 
in either direction from this estimate cannot be considered as “incorrect.” Unfortunately, tests with different 
time steps rarely give any strong indications toward an optimal value. 7 Thus, some ambiguity remains in 
the selection of the time step size. Solutions were therefore obtained on the baseline mesh using the smaller 
time step to test the effects of reducing the time step size for a fixed mesh. 

E. Unstructured Mesh Generation 

Three different meshes were generated for the extruded wing configuration: a relatively coarse mesh that 
was used solely for RANS computations, a baseline mesh, and a refined mesh. The coarse RANS mesh 
was generated using SolidMesh. SolidMesh is an interface to the unstructured surface and volume mesh 
generation software AFLR2 and AFLR3. 25 AFLR3 uses an advancing front algorithm to insert a point in 
the mesh. The point insertion is followed by a local reconnection to improve mesh quality. The baseline 
and refined meshes were generated using Grid Tool 2 6 and VGridns. 27 VGridns uses an advancing layer 
algorithm to generate a tetrahedral volume mesh. The near-body elements of these all-tetrahedral meshes 
were converted into prisms using the Blacksmith utility 28 thereby producing a mixed element hybrid mesh. 
SolidMesh produces a hybrid mesh with prisms and tetrahedra automatically. These hybrid meshes were 
employed because of their potential for improved efficiency and accuracy in comparison to unstructured 
tetrahedral meshes. The GridTool/ VGridns combination was chosen for the DES meshes because it gives 
the user control of the point spacing on the wing surface and the flow field through the use of “sources.” 

Figure 2 shows two cross-sections through the refined DES mesh. Since the figure was generated using 
cutting planes, the line segments represent intersections of cell faces with the cutting plane. The connec- 
tivity of the mesh is evident and the tetrahedra/prism layers near the surface are clearly visible. Note the 
finer resolution in the focus region downstream of the horn for the refined mesh. This figure also illustrates 
the connection between the surface mesh and volume mesh. Because of the manner in which the mesh is 
generated, i.e., anisotropic tetrahedra/prisms transitioning to isotropic tetrahedra, the surface mesh charac- 
teristics are propagated into the volume mesh. Therefore, in order to refine the mesh spacing in the domain, 
it is necessary to refine the mesh on the surface. 

Mesh refinement was obtained by reducing the magnitudes of the line sources in appropriate regions. 
In this case, the mesh was locally refined by reducing the source strengths in the upper surface FR by a 
factor of two while holding the other source strengths constant. This was found to be a reasonably efficient 
mechanism for mesh refinement for this problem. Since the focus of this effort was the flow on the upper 
surface, the lower surface mesh was made relatively coarse in comparison. 

Table 1 shows statistics for the three meshes employed to generate the results reported here. Notice 
that the localized mesh refinement in the refined mesh produces more than a 25% increase in the number 
of cells when compared to the baseline mesh. In some sense, this mesh refinement can be thought of as a 
manual, albeit crude, attempt at solution adaptive meshing. In all cases, the distance to the first point off 
the wall was defined so that an average y + value less than 0.5 was obtained. This value is well within the 
recommended values for the turbulence model employed. 
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Figure 2. Refined DES mesh in cutting planes for extruded GLC305/944-ice shape (extruded wing) 


Table 1. Statistics for the meshes employed in this effort 



# Nodes 

# Faces 

# Cells 

Extruded Wing (coarse mesh) 

890,958 

5,299,738 

2,258,447 

Extruded Wing (baseline mesh) 

1,828,711 

17,327,367 

8,341,019 

Extruded Wing (refined mesh) 

2,394,393 

22,088,508 

10,579,834 


F. Flow Solution 

The flow solver employed in this effort is the non-commercial version of Cobalt 60) 29 which is now called AVUS 
to reduce confusion with the commercial version of the flow solver. The AVUS flow solver was designed for 
general unstructured meshes. It employs a nonlinear Riemann solver for the inviscid flux computations and 
can be run either in explicit or implicit mode. Second-order spatial accuracy is obtained using a linear least- 
squares reconstruction of the data. The detailed computational methodology used in AVUS is described in 
the literature. 18 

All computations were performed using second-order spatial accuracy. The RANS solutions reported here 
were obtained using first-order, implicit local time stepping and do not in any way represent time-accurate 
solutions. For the DES computations, the second-order temporal integration scheme was employed. It should 
be noted that a CFL number of unity in the focus region yields a CFL number on the order of 100,000 in 
the boundary layer. Therefore, for stability reasons, it was necessary to use the fully-implicit integration 
scheme (0—1) for the DES computations. Based on recommendations in the Cobalt 6 o user’s manual, 28 
two Newton iterations, each consisting of 30 Gauss-Seidel iterations, were employed per time step. Several 
turbulence models are available including the Spalart-AIlmaras one-equation model, 19 which was used for 
the computations reported here. No transition point was specified, and the flow was assumed to be fully 
turbulent. A slip boundary condition was applied on the artificial side boundaries to reduce computational 
expense. This boundary condition is not believed to introduce significant error in the modeling of the flow 
since a porous sidewall was employed for boundary layer control in the experiment. 


V. Results 

In this section, we present comparisons between experimental data 16 ’ 17 and results computed using both 
steady RANS and unsteady DES models. We focus first on gross flow field quantities, e.g., lift and drag 
coefficients. We then present comparisons of detailed flow field quantities at the midspan of the wing - the 
mean surface pressure distribution and field distributions of the mean streamwise velocity component and 
the RMS of the streamwise velocity fluctuations. We conclude with an assessment of the validity of the DES. 
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More details and additional comparisons are provided in Mogili 30 and Thompson and Mogili. 31 

Because of the associated computational expense, we have chosen to apply DES only to the 4 deg and 
6 deg angle of attack cases. RANS computations were performed for angles of attack from 0 deg to 6 deg in 
2 deg increments. Although the lower surface also has a region of separated flow, the flow is dominated by 
the large region of unsteady separated flow on the upper surface. Additionally, no field data were obtained 
on the lower surface during the experiments. Therefore, our discussion will focus on the region of the flow 
downstream of the upper surface horn. 

The RANS computations were performed on three meshes - the coarse mesh, the baseline mesh, and the 
refined mesh - described in Table 1. The DES computations were performed on the baseline mesh and the 
refined mesh. The RANS simulations were initiated using an impulsive start and were continued until the 
normal force reached a steady state, which typically required 15,000-16,000 global iterations. The RANS 
solutions were then used as initial conditions for the DES computations. The sudden change in the eddy 
viscosity (due to the reduction in the length scale when the DES computation is initiated) represents a non- 
physical transient that must be eliminated. Therefore, in all cases, DES computations were performed for 
0.75s (a nondimensional time of 17 00 t/L=33.7), which corresponds to 5000 steps for Ai=0.15ms and 10000 
steps for 075ms to eliminate this transient. Data was then collected during an additional 0.75s for a 
total nondimensional time of 67.4 beyond the steady RANS solution. 

All cases reported in this effort were run on 64 processors on the EMPIRE cluster or the MAVERICK 
cluster at the ERC at Mississippi State University. The EMPIRE cluster is a supercomputer class cluster 
of workstations consisting of 1038 one GHz or better Pentium III processors each with one or more GB of 
RAM. The MAVERICK cluster is also a supercomputer class cluster comprised of 192 IBM x335 nodes. 
Each node contains dual 3.06 GHz Xeon processors and 2.5 GB of RAM. For the three-dimensional RANS 
solutions on the baseline mesh, 1 hour and 39 minutes were required for each 100 iterations on the EMPIRE 
cluster, whereas 2 hours and 23 minutes were required for each 100 iterations on the refined mesh. For 
the three-dimensional DES computations, 2 hours and 52 minutes were required for 100 time steps on the 
baseline mesh. Three hours and 52 minutes were required for each 100 time steps on the refined mesh. 

A. Description of Experimental Data 

Addy, et al. 16 performed icing effects studies for the configuration of interest in this effort. The lift and 
pitching moment data were obtained by integrating the surface pressures, while the drag coefficients were 
calculated using the standard momentum deficit method based on pressures measured using a wake probe. 
Corrections to the integrated performance coefficients accounting for solid and wake blockage and streamline 
curvature were applied to the data during post processing. Broeren, et al 17 also carried out flow field 
measurements on the upper surface of the same model. Data were obtained at three different angles of 
attack preceding stall at Reynolds numbers of 3.5xl0 6 and 6.0xl0 6 and Mach numbers of 0.12 and 0.21. 
Spht-hot-film anemometry was used to measure the time-averaged flow velocities and its RMS fluctuations. 

B. Comparison of Lift and Drag Coefficients 

Figure 3 shows a comparison of the predicted lift and drag coefficients with experimental data. For the DES 
results, the lift and drag coefficients represent mean values that were obtained using the temporal averaging 
procedure described above. In general, the steady RANS methodology underpredicts the lift with respect to 
the experimental data with the differences increasing with increasing angle of attack. Mesh refinement does 
not significantly improve the RANS results. The DES results show significant improvement in comparison 
to e steady RANS results at an angle of attack of 6 deg, irrespective of the mesh density or time step 
size Less improvement is apparent at a 4 deg angle of attack. With the exception of the steady RANS 

° b ! amed ° n the coarse mesh > the dra g coefficient comparison shows a similar trend. The steady 
RANS and time-averaged DES approaches underpredict the drag coefficient although the DES results show 
improvements relative to the RANS results. The relatively good agreement between the experimental data 
and the steady RANS results computed on the coarse mesh probably should be considered fortuitous since 
the agreement degrades as the mesh is refined. 
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Angle of attack (deg) 


(a) Lift coefficients 

Figure 3. Comparison of predicted wing lift 



Angle of attack (deg) 


(b) Drag coefficients 

drag coefficients and experimental data 16, 17 


C. Detailed Comparisons with Experimental Data 

We now present comparisons of detailed flow field data for the predicted results and experimental data at 
angles of attack of 4 deg and 6 deg. In general, there were significant differences between steady RANS 

results and time-averaged DES results. The effects of spatial and temporal refinement were less significant 
for a given method. 

Figure 4 shows a comparison of predicted results and experimentally determined time-averaged chord- 
wise pressure coefficient distributions at the midspan of the wing. Figure 4(a) shows comparisons between 
experimental data, steady RANS results, and time-averaged DES results for an angle of attack of 4 deg. 
In general, none of the predicted results do a very good job of reproducing the upper surface pressures 
downstream of the horn. The constant pressure plateau appearing just downstream of the horn in the DES 
results is indicative of a significant region of flow separation. The pressure plateau is less prominent in 
the experimental data and is not evident in the RANS results. The pressure recovery in the DES results 
shows an aft shift in comparison to the RANS results and better matches the experimental data. We note 
that, with the exception of the coarse RANS solution, spatial and temporal refinement have only minimal 
impact on the resulting pressure distributions. The lift values obtained by integrating the pressures appear 
consistent with the results observed in Figure 3. 

Figure 4(b) shows the same comparisons for the 6 deg angle of attack case. Again, we note that the 
discrepancies between the steady RANS prediction made on the coarse mesh and the experimental data 
are decreased as the mesh is refined. Less influence of spatial and temporal refinement is observed for 
o er results. The pressure plateau exhibited in each distribution suggests an extensive region of separated 
flow. The pressure recovery displayed in the time-averaged DES results shows better agreement with the 
experimental data than the results obtained from the steady RANS simulations. As evident from the plot 
ere is a downward ’ shift of the predicted DES pressures relative to the experimental data. The shape 
of the pressure distribution is well predicted except on the lower surface just aft of the horn for which 
the resolution may be inadequate since we chose to focus on the upper surface, and just downstream of 
the horn. This is particularly true for the refined mesh case. As we will note below, there is a significant 
secondary recirculation in the predicted results whose streamwise extent appears to coincide with the location 
of the differences between the predictions and the experimental data on the upper surface of the wing The 
qualitative agreement between the overall shape of the time-averaged DES pressure distributions and the 

experimental data is encouraging as it represents an improvement when compared to the steady RANS 
predictions. J 
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Figure 4. Comparison of predicted midspan wing pressure coefficients and experimental data 


We now consider the distribution of the mean streamwise velocity in a plane parallel to the chordline 
at the midspan of the wing. Due to space limitations, we focus on results obtained on the refined mesh. 
However, the comparisons shown here are representative of those observed during the study. Figure 5(a) 
shows a comparison between the predicted midspan streamwise velocity contours and experimental data 
obtained using a hot split-film probe for the refined mesh for an angle of attack of 4 deg. The location of the 
reattachment of the primary upper surface flow at this spanwise position may be estimated by finding the 
position at which the zero streamwise velocity contour intersects the surface of the wing. From Figure 5(a) , we 
on£ Se ! t u at r i eattachment occurs at approximately 31% of the chord in the RANS results, at approximately 
39% of the chord for the DES results, and approximately 28% of the chord for the experimental results. 
Thus, for both RANS and DES results, the predicted reattachment location is aft of the location suggested 
by the experimental data. This is a trend we observed for all of the computed results regardless of the 
methodology employed or spatial and temporal resolution. Further, the DES results show a rearward shift of 
the region of maximum reversed flow velocity in comparison to both the RANS results and the experimental 
data. Additionally there is a significant secondary recirculation evident in the DES results just downstream 
0t *7 e (as indicated by the green region near the surface). There is a very small recirculation region 
m the RANS results that is not visible in Figure 5(a). Although it is likely that a recirculation of lesser 
extent exists m the flow field, it does not appear in the experimental data. The seemingly anomalous velocity 
contours that appear m the experimental data just downstream of the horn are artifacts from the process 
employed to generate the contour plots. 

Figure 5(b) shows a comparison between the predicted midspan streamwise velocity contours and experi- 
mental data obtained using a hot split-film probe for the refined mesh for an angle of attack of 6 deg. For this 
case, both the RANS and DES results show a fully-separated flow on the upper surface. The experimental 
data suggests a reattachment location at approximately 50% chord. There is significantly more reversed 
° W no 6 7 ^ NS resuits in comparison to the DES results as evidenced by the thick contour of velocity in 
the -0.2 to 0.0 range that extends to the trailing edge in the RANS results. The RANS results also show a 
reversed flow velocity of larger magnitude from approximately 5% to 50% of the chord. As before, the region 
of maximum reversed flow is shifted aft in the predicted results relative to the experimental, more so in the 
DES results than in the RANS results. Again, a significant secondary recirculation is present in the DES 
results just downstream of the horn. Although not evident from this image, there is a smaller recirculation 
region present m the steady RANS results. At this angle of attack, there is a small recirculation region is 
barely perceptible in the experimental data. 

We now consider the streamwise velocity fluctuations in the midspan plane. Figure 6(a) shows a com- 
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Time-averaged DES results (ZU=0.075ms) Time-averaged DES results (Af=0.075ms) 



Experimental hot-split-film data 17 Experimental hot-split-film data 17 

(a) 4 deg angle of attack (b) 6 deg angle of attack 


Figure 5. Comparison of predicted midspan u-velocity contours and experimental data for the extruded wing 
— refined mesh 


parison between the RMS of the fluctuations in the streamwise velocity component predicted by the DES 
computation and experimental data. The turbulence intensity was calculated as the root-mean-square of 
the fluctuating streamwise component of velocity normalized by the freestream velocity. The fields are 
qualitatively similar. However, the location of the region containing the maximum fluctuations is shifted 
downstream. Additionally, there is a region in the predicted results just downstream of the horn that is 
characterized by a much lower level of fluctuations. The secondary circulation mentioned above is located 
in this region. These results suggests that the shear layer roll up may be delayed. 

Figure 6(b) shows a comparison between the RMS of the fluctuations in the streamwise velocity com- 
ponent predicted by the DES computation and experimental data. The fields are qualitatively similar. 
However, the location of the region containing the maximum fluctuations is shifted downstream and its 
downstream extent is somewhat larger than that observed in the experimental data. This suggests that the 
predicted upper surface flow is somewhat more active than the physical flow. As before, the region of lower 
values of predicted fluctuations in the region just downstream of the horn suggests that the previously noted 
recirculation is a fairly steady phenomenon. Again, the shear layer roll up appears somewhat delayed in 
comparison to the experimental data. 

D. Validity of DES Results 

It is now useful to question the validity of the DES. In this context, validity refers to how well the DES has 
resolved significant features in the flow, i.e, are the spatial and temporal resolutions adequate for the problem 
at hand? The approach taken here is to compare solutions computed on two different meshes - the baseline 
and refined meshes which differ by the spacing of mesh points in the focus region — and solutions computed 
on the same mesh using two different time steps - the refined mesh with At=0.15ms and At=0.075ms. 
Because of the complexity of the unsteady three-dimensional flow fields, the comparisons employed here are 
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DES results (zlt=0.075ms) 


DES results (.4£=Q. 075ms) 





Experimental hot-split-film data 17 
(a) 4 deg angle of attack 


Figure 6. Comparison of the RMS of the fluctuate 
refined mesh 


^/tr 



Experimental hot-split- film data 17 
(b) 6 deg angle of attack 

in the u-velocity component for the extruded wing — 


based on the lift histories. The periodogram function of MATLAB 32 was used to analyze the time histories 
of the wing lift coefficient. The power spectral density (PSD) of the time history of the lift coefficient is 
computed and plotted against the frequency of the signal. Due to space limitations, we focus only on the 
results obtained after mesh refinement. The time interval was 0.75s to 1.5s after the initiation of the DES 
computation. For ZU=0.075ms, this corresponds to 10,000 time steps. For ^t=0.15ms, this corresponds 
to 5,000 time steps. Figure 7 shows a comparison of the PSD signatures for the three-dimensional DES 
solutions on the baseline and refined meshes with Zlf=0.075ms for 6 deg angle of attack. From the figure, 
it is apparent that similar signatures are obtained up to 15Hz. Figure 7(b), which is the zoomed version of 
Figure 7(a) in the frequency range 10-100Hz, shows differences above 15Hz. Similar results were obtained for 
the time step refinement study but with quantitative differences occurring at frequencies of approximately 
12Hz. 

These and similar results indicate that both the time-step refinement and mesh refinement only signif- 
icantly affected the signal at frequencies higher than 12-15Hz. The power at these frequencies is at least 
an order of magnitude less (range 10 04 to 10 -06 ) than the power at the lower frequencies. Taken together 
with the fact that the time-averaged quantities are so similar, we can conclude that these refinements did 
not significantly affect the results of the simulations. Since the character of the solutions did not change 
with these refinements, this suggests that we have valid DES results within the context of the AVUS flow 
solver. This does not address the accuracy of the simulation. 

We also investigated the effects of changing the number of Newton iterations employed for each time 
step. We performed an AVUS computation for the 6 deg angle of attack condition on the refined mesh with 
a time step of 0.075ms. We employed three Newton iterations instead of two. The results obtained using 
three Newton iterations per time step are quantitatively very similar to those obtained using two iterations. 
In fact, the differences were within those obtained after spatial or temporal refinement. We attribute this 
to the fact that we have selected a time scale that is appropriate for the problem under consideration. We 
also performed a PSD analysis similar to those described above. Again, we found the differences between 
the three iteration and two iteration cases to be similar to the differences that were observed for spatial and 
temporal refinement. 


VI. Conclusion 

The problem of flow field simulation for iced wing configurations is a complex one that severely taxes 
existing capabilities for geometry modeling, mesh generation, and flow solution. Particularly challenging are 
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Figure 7. Comparison of power spectral density plots 
attack, with Zl £=0.075 ms 



for baseline mesh and refined mesh at 6 deg angle of 


cases in which there are regions of unsteady flow separation. In this context, the effectiveness of Detached- 
bddy Simulation (DES) as a tool for predicting icing effects was evaluated. The AVUS code was employed 
iDA°^ Ute solutlons for an iced win S configuration using DES and steady Reynolds Averaged Navier-Stokes 
(KANa) equation methodologies. The configuration was an extruded GLC305/944-ice shape section with a 
rectangular planform. The model was mounted between two walls so no tip effects were considered. The 
numerical results were validated by comparison with experimental data for the same configuration. 

h e benefits of employing three-dimensional DES computations are not clear at this time. The time- 
averaged DES computations did show improvement in lift and drag results at an angle of attack of 6 dex 
when compared to steady RANS results. Similar agreement was not obtained at an angle of attack of 4 deg. 
Comparisons of detailed flow field quantities such as surface pressure distribution and the mean streamwise 
flow velocity distribution did not show the level of agreement suggested by the integrated quantities. In 
particular DES results showed more extensive flow separation than the experimental data. Further, the 
ES results showed a significant secondary recirculation that was not present in the experimental data. The 
D.ho results also exhibited a delayed shear layer roll up. 

Based on our results, we believe that DES may prove useful in a limited sense to provide analysis 
of iced wing configurations when there is significant flow separation, e.g., near stall, where steady RANS 
computations are demonstrably ineffective. However, more validation is needed to determine what role DES 
can play as part of an overall icing effects prediction strategy. 


VII. Acknowledgment 

This effort was supported by NASA GRC (NAG3-2562 and NAG3-2892) with Yung Choo as Technical 
Monitor We would like to extend our appreciation to Stuart Pope of NASA LRC for his assistance with 
GndTool and Matthew Grismer of AFRL for his assistance with the AVUS flow solver. We would also like to 
h an k An dy Broeren of UIUC for his helpful discussions relating to the experimental data. Thanks are also 
due Mithun Varma of MSU who assisted with the images and Satish Chalasani who generated the coarse 


References 


^Bragg M B., Broeren, A. P., and Blumenthal, L. A., “Iced-Airfoil and Wing Aerodynamics,” SAE, 
3 * ani > *‘5 Low Speed Flows Involving Bubble Separations, Pergamon, New York, 1964. 

87-010i° taPCZUk ’ M " " NumeriCaI Analysis of aa NACA0012 Airfoil with Leading Edge Ice Accretions,” 


June 2003. 

1987, AIAA Paper 


12 of 13 


American Institute of Aeronautics and Astronautics 





of im P A^LZ d ^0sk. “ Navi0r - Stokes Analysis of the Flowed Characteristics 

AIAA Paper 1999-0093. ’ ’ EffeCtS ° f S,mula ted-Spanwise-Ice Shapes on Airfoils: Computational Investigation,” 1999, 

rspWrt’ P°r.' -P agg ’ M n B '.; ™ SimuIations ^oOm with Ice Shapes,” 2003, AIAA Paper 2003-0729. 

July 2001. ’ ’ ^ erSOnS m e to etached-Eddy Simulation Grids,” Tech. Rep. NASA CR-2001-211032, NASA, 

Hybrid ffi/lBS Ap^roIIh,” ^DNS/LEsTt AFOSRfruZT Z o' ° f LES f ° r Wn8s ’ “ d °” a 

Forsythe, J. R„ Squires K D Wurtt O / \ Conference on DNS/LES, August 1997. 

High Alpha,” 2002, AIAA Paper 2002-0591. ’ ’ ’’ Spalart > p - R-> Detached-Eddy Simulation of Fighter Aircraft at 

Eddy Emulation H%^^ynoids^umber fep^rated Hoira^^rMeedi E " “fTZZT ° f Unstructured Grids for Detached 
Gnd Generation in Computational Field Simulations, Honolulu, June 2002. O/ ‘ Internatlonal Conference on Numerical 

J - 

t Kumar, S. and Loth, E "Detached Eddy Simulations of an Iced-Airfoil,” 2001 AIAA Pant 20oTo678 

-Tromps:: D Moth? ta ctlt,fs S 11;tH f ° r d Ch" Wi Y ke Shapes> ” 2C04 ’ AIAA Papw 2004 - 0564 - 

4dt: T X NASA 

Paper 2003-0727. ’ ’’ d Tunnel Stud y of Icing Effects on a Business Jet Airfoil,” 2003, AIAA 

2004! MtX Papt tot y 0 5 H 59 E '' ^ ®- ““' d M *— ab °^ an Airfoil with Leading-Edge Ice Shapes,” 

Euler/Navier-Stokes Flow Solver," ’ IMS,' AllTpa^er i99?0786 fimng Meth ° ds ° f Cobalt6 °- A Parallel, Implicit, Unstructured 
tiale, Voh'l S ' R '' " A 0ne - Ec i uation Turbulence Model for Aerodynamic Flows,” La Recherche Aerospa- 

with^Leading-Edge Ice Accretion,” 2004, AIAA Papertooio'sffl' ° FD AnalySIS ° f the Aer °dynamics of a Business Jet Airfoil 
Flow with an Unstructured sXer^ 2000^1^ Eddy SimuIation of a Supersonic Axisymmetric Base 

WeatheTRev'^l 91,’ lGesTp^MM^™ Experim6ntS with the Primitive Equations. I. The Basic Experiments,” Man. 

f 7 ICe , A ; Cr f° n Predictio “.” 2 °03, AIAA Paper 2003-1070. 

put “j^ m H l :: b zc°g^ s ofthe 6th in Com - 

Grid -“T^t“t pn °t “ d CRcVrts,l» ^VttnlTd U “ red 

elias 4t“ 

No. 1, January 1996, pp. 43-^“ enS10Ila nStrUCtUred Vlscous Grlds b X Advancing Layers Method,” AIAA Journal, Vol. 34, 
“Strang, W. Z., Cobalt 6Q : User’s Manual , September 2000. 

structured Meshes,” 1997, AU^ Paper^^asS. L ^ ' ImpllClt AIg ° rithm for Solvin S Time Dependent Flows on Un- 
AugSo4.' P ” RANS and DES CWp “^'™ f0r Wm ° S Accretions, Master’s thesis, Mississippi State University, 

Report (NAG3-2892),” tt^rep ^NA^ Attttetemtr^OOA ° f Separated FIow ^nd Wings with Ice Accretions: Year 1 
2 http :/ /www . math works . com/ . 


13 of 13 


American Institute of Aeronautics and Astronautics 



